Hazard Modeller's Toolkit - Basic Catalogue Tools

The following is a demonstration of the basic ancilliary functions of the HMTK for the Catalogue Workflow

The demonstration will do the following:

1) Load the Earthquake Catalogue

2) Explore Basic Functions for Investigating and Visualising the Properties of the Catalogue

(Python Step) Import the needed libraries



In [1]:

    
# Python Numerical and Plotting Libraries
import numpy as np
import matplotlib.pyplot as plt

# HMTK Catalogue Import/Export Libraries
from hmtk.parsers.catalogue.csv_catalogue_parser import CsvCatalogueParser, CsvCatalogueWriter

# HMTK Plotting Tools
from hmtk.plotting.seismicity.catalogue_plots import (plot_depth_histogram,
                                                      plot_magnitude_time_scatter,
                                                      plot_magnitude_time_density,
                                                      plot_magnitude_depth_density,
                                                      plot_observed_recurrence)
from hmtk.plotting.mapping import HMTKBaseMap
print 'Imports OK!'









    



Imports OK!

Import the Catalogue



In [2]:

    
input_catalogue_file = 'input_data/Aegean_ExtendedCat1.csv'
parser = CsvCatalogueParser(input_catalogue_file)
catalogue = parser.read_file()
print 'Input complete: %s events in catalogue' % catalogue.get_number_events()
print 'Catalogue Covers the Period: %s to %s' % (catalogue.start_year, catalogue.end_year)









    



Catalogue Attribute Identifier is not a recognised catalogue key
Input complete: 16401 events in catalogue
Catalogue Covers the Period: 1904 to 2012

If the catalogue is not in chronological order this can cause errors in later functions - sort chronologically after import to fix this



In [3]:

    
# Sort catalogue chronologically
catalogue.sort_catalogue_chronologically()
print 'Catalogue sorted chronologically!'









    



Catalogue sorted chronologically!

Viewing the Catalogue



In [4]:

    
# Configure the limits of the map and the coastline resolution
map_config = {'min_lon': 18.0, 'max_lon': 32.0, 'min_lat': 33.0, 'max_lat': 43.0, 'resolution':'h'}
# Create a hmtk basemap
basemap1 = HMTKBaseMap(map_config, 'Earthquake Catalogue')
# Add a catalogue
basemap1.add_catalogue(catalogue)

Adding a Source Model

The following is an area source model for the Aegean region (Weatherill & Burton, 2010) - this source model will be used in further examples



In [5]:

    
# Source Model file
source_model_file = 'input_data/WT2006_Aegean_Sources.xml'

# Import the source model parser
from hmtk.parsers.source_model.nrml04_parser import nrmlSourceModelParser

# Load source model
parser = nrmlSourceModelParser(source_model_file)
source_model = parser.read_file()

# Create a hmtk basemap with a catalogue and source model
basemap2 = HMTKBaseMap(map_config, 'Earthquake Catalogue & Area Source Model')
basemap2.add_catalogue(catalogue, overlay=True) # Adds the catalogue
basemap2.add_source_model(source_model, area_border='r-', border_width=1.5) # Adds the source model









    



Area source - ID: 1, name: Montenegro
Area source - ID: 2, name: Albania
Area source - ID: 3, name: TRGZ
Area source - ID: 4, name: Cephalonia Transform
Area source - ID: 5, name: South Ionian
Area source - ID: 6, name: Gulf of Corinth
Area source - ID: 7A, name: West Hellenic Outer
Area source - ID: 7B, name: East Hellenic Outer
Area source - ID: 8A, name: West Hellenic Inner
Area source - ID: 8B, name: East Hellenic Inner
Area source - ID: 9, name: Cyclades
Area source - ID: 10, name: North Anatolian Fault - Main Branch
Area source - ID: 11, name: North Anatolian Fault - Southenrn Branch

Magnitude-Time Properties



In [6]:

    
# Show distribution of magnitudes with time
plot_magnitude_time_scatter(catalogue, fmt_string='.')



In [7]:

    
# Add some error-bars
plot_magnitude_time_scatter(catalogue, plot_error=True, fmt_string='.')



In [8]:

    
# Plot the magnitude-time density
magnitude_bin = 0.1
time_bin = 1.0 # in Decimal Years
plot_magnitude_time_density(catalogue, magnitude_bin, time_bin)

Simple Histograms - Depth



In [9]:

    
# Depth histogram
plot_depth_histogram(catalogue, 5.)



In [10]:

    
# Depth histogram (with bootstrap sampling depth uncertainty)
depth_bin = 5.0 # in km
plot_depth_histogram(catalogue, depth_bin, bootstrap=500)









    



/home/gw/Documents/hmtk/hmtk/hmtk/seismicity/utils.py:275: RuntimeWarning: divide by zero encountered in divide
  lower_bound = (bounds[0] - data) / uncertainties



In [11]:

    
# Magnitude Depth Density
magnitude_bin = 0.1
depth_bin = 5.0  # in km
plot_magnitude_depth_density(catalogue, magnitude_bin, depth_bin, logscale=True)

Simple plots of earthquake recurrence



In [12]:

    
# Time-varying completeness
completeness = np.array([[1985., 4.0],
                         [1964., 5.0],
                         [1910., 6.5]])
plot_observed_recurrence(catalogue, completeness, 0.1, catalogue.end_year)

Using Numpy logical tools for simple selection process



In [13]:

    
# Limit the catalogue to the time period 1960 - 2012
valid_time = np.logical_and(catalogue.data['year'] >= 1960,
                            catalogue.data['year'] <= 2012)
catalogue.select_catalogue_events(valid_time)
plot_magnitude_time_density(catalogue, 0.1, 1.0)
print 'Catalogue now contains %s events' % catalogue.get_number_events()









    












    



Catalogue now contains 16282 events



In [14]:

    
# Limit the catalogue to depths less than 50 km
valid_depth = catalogue.data['depth'] <= 50.
catalogue.select_catalogue_events(valid_depth)
plot_depth_histogram(catalogue, 2.0)

Exporting the Catalogue

At any time the current content of the catalogue can be written to csv:



In [15]:

    
# Set-up the file writer
output_file_name = 'output_data/basic_demo_catalogue_1.csv'
writer = CsvCatalogueWriter(output_file_name)

# Write the catalogue to file
writer.write_file(catalogue)

print 'File %s written' % output_file_name









    



File output_data/basic_demo_catalogue_1.csv written

If you wish, you can export only the complete events:



In [16]:

    
completeness = np.array([[1985., 4.0],
                         [1964., 5.0],
                         [1910., 6.5]])
# Set-up the exporter
output_file_name = 'output_data/basic_demo_catalogue_complete_1.csv'
writer = CsvCatalogueWriter(output_file_name)

# Write the catalogue to file, purging events from the incomplete period
writer.write_file(catalogue, magnitude_table=completeness)

print 'File %s written' % output_file_name









    



[False False False ..., False False False]
[False False False ..., False False False]
[False False False ..., False False False]
File output_data/basic_demo_catalogue_complete_1.csv written